# Analyzing metabolomics data for Mutographs ----
# 04_random_forests
# Running random forest models

rm(list = ls())

# 0 - Definition of libraries, paths and functions ----

Sys.setlocale("LC_TIME", "C")

library(stringr)
library(tidyr)
library(Hmisc)
library(MASS)
library(data.table)
library(dplyr)
library(forcats)
library(tibble)
library(openxlsx)
library(ggplot2)
library(cowplot)
library(scales)
library(patchwork)
library(ggpubr)

# Loading previously processed data
metabo_PROC <- readRDS("output/metabo_PROC")
metabo_MAN <- readRDS("output/metabo_MAN")
metabolites <- names(metabo_PROC)[str_detect(names(metabo_PROC),"@")]

cosmic_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_SBS96_abs_mutations.csv") %>% 
  rename(donor_id=X, SBS1536A_cosmic = SBS1536A, SBS1536B_cosmic = SBS1536B, SBS1536F_cosmic = SBS1536F, SBS1536I_cosmic = SBS1536I)

# Adding mutational burden
mutburden_sbs <- read.csv("data/sigs/output_RCC_Manuscript_COSMIC_SBS96_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden_sbs=Mutational.burden) %>% select(donor_id, mutburden_sbs)
mutburden_dbs <- read.csv("data/sigs/output_RCC_Manuscript_COSMIC_DBS78_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden_dbs=Mutational.burden) %>% select(donor_id, mutburden_dbs)
mutburden_id <- read.csv("data/sigs/output_RCC_Manuscript_COSMIC_ID83_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden_id=Mutational.burden) %>% select(donor_id, mutburden_id)
mutburden_cn <- read.csv("data/sigs/output_RCC_Manuscript_COSMIC_CNV48_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden_cn=Mutational.burden) %>% select(donor_id, mutburden_cn)
mutburden_sv <- read.csv("data/sigs/output_RCC_Manuscript_denovo_SV32_stat_info.csv") %>% 
  rename(donor_id=Sample, mutburden_sv=Mutational.burden) %>% select(donor_id, mutburden_sv)
mutburden <- mutburden_sbs %>% left_join(mutburden_dbs) %>% left_join(mutburden_id) %>% left_join(mutburden_cn) %>% left_join(mutburden_sv)

mutburden_all <- names(mutburden)[2:6]

metabo_MAN_burden <- metabo_MAN %>% left_join(mutburden)

# Checking AA and SBS12 exposed cases
# Listing cases with >10% of signature attribution due to SBS22a + SBS22b or SBS12
cosmic_relative <- cosmic_attrib %>% 
  mutate(total_attrib = rowSums(.[2:15])) %>% 
  mutate(across(starts_with("SBS"), ~./total_attrib))

aa_cases <- filter(cosmic_relative, SBS22 + SBS1536I_cosmic > 0.1)$donor_id
SBS12_cases <- filter(cosmic_relative, SBS12 > 0.1)$donor_id

metabo_MAN_burden_rf <- metabo_MAN_burden %>% filter(!donor_id %in% c(aa_cases,SBS12_cases))

# For SVs de novo
sv_dn_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_denovo_SV32_abs_mutations.csv") %>% rename(donor_id=X)
sv_dn <- names(sv_dn_attrib)[2:4]
metabo_MAN_SV <- metabo_MAN %>% left_join(sv_dn_attrib) %>% # Adding SV DN signatures
  mutate(across(
    any_of(c(sv_dn)),
    .fns = list(cat = ~cut2(.,g=2), int = ~qnorm((rank(.,na.last="keep")-0.5)/sum(!is.na(.))), logdelta = ~log2(.+1)),
    .names = "{.col}_{.fn}"))

# Running random forest models for interesting signatures ----
# Same as regressions: if sparseness of signature is > 30%, a classification model is ran
# Else, a regression model is chosen

# In all countries
library(randomForest)

set.seed(42)
# SBS4 - Categorical
rf_SBS4 <- randomForest(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                        y = metabo_MAN$SBS4_cat)
randomForest::varImpPlot(rf_SBS4)

# SBS22a - Categorical
rf_SBS22a <- randomForest(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                          y = metabo_MAN$SBS22_cat)
randomForest::varImpPlot(rf_SBS22a)

# SBS22b - Categorical
rf_SBS22b <- randomForest(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                          y = metabo_MAN$SBS1536I_cosmic_cat)
randomForest::varImpPlot(rf_SBS22b)

# SBS40b - Continuous
rf_SBS40b <- randomForest(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                          y = metabo_MAN$SBS1536A_cosmic)
randomForest::varImpPlot(rf_SBS40b)

# DBS2 - Categorical
rf_DB2 <- randomForest(x = select(filter(metabo_MAN, !is.na(DBS2)), c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                        y = filter(metabo_MAN, !is.na(DBS2))$DBS2_cat)
randomForest::varImpPlot(rf_DB2)

# ID1 - Continuous
rf_ID1 <- randomForest(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                          y = metabo_MAN$ID1)
randomForest::varImpPlot(rf_ID1)

# ID5 - Continuous
rf_ID5 <- randomForest(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                          y = metabo_MAN$ID5)
randomForest::varImpPlot(rf_ID5)

# ID8 - Continuous
rf_ID8 <- randomForest(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                          y = metabo_MAN$ID8)
randomForest::varImpPlot(rf_ID8)

# Only in Romania and Serbia for SBS22-like signatures
# SBS22a - Categorical
rf_SBS22a_sro <- randomForest(x = select(filter(metabo_MAN, country %in% c("Romania","Serbia")), c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                              y = filter(metabo_MAN, country %in% c("Romania","Serbia"))$SBS22_cat)
randomForest::varImpPlot(rf_SBS22a_sro)

# SBS22b - Categorical
rf_SBS22b_sro <- randomForest(x = select(filter(metabo_MAN, country %in% c("Romania","Serbia")), c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                              y = filter(metabo_MAN, country %in% c("Romania","Serbia"))$SBS1536I_cosmic_cat)
randomForest::varImpPlot(rf_SBS22b_sro)

# Saving all RF models in a single list
l_rf <- list()
l_rf[["SBS4"]] <- rf_SBS4
l_rf[["SBS22a"]] <- rf_SBS22a
l_rf[["SBS22b"]] <- rf_SBS22b
l_rf[["SBS40b"]] <- rf_SBS40b
l_rf[["SBS22a_sro"]] <- rf_SBS22a_sro
l_rf[["SBS22b_sro"]] <- rf_SBS22b_sro
l_rf[["DBS2"]] <- rf_DB2
l_rf[["ID1"]] <- rf_ID1
l_rf[["ID5"]] <- rf_ID5
l_rf[["ID8"]] <- rf_ID8
# saveRDS(l_rf, "results/random_forests/rf_all")

# Running RF permutation to calculate p-values for each features
library(rfPermute)

set.seed(42)
# With 1000 permutations for each model of interest
rp_SBS4 <- rfPermute(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                     y = metabo_MAN$SBS4_cat, # Categorical
                     num.rep=1000)
# saveRDS(rp_SBS4, "results/random_forests/rp_SBS4")

rp_SBS13 <- rfPermute(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                       y = metabo_MAN$SBS13_cat, # Categorical
                       num.rep=1000)
# saveRDS(rp_SBS13, "results/random_forests/rp_SBS13")

rp_SBS18 <- rfPermute(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                      y = metabo_MAN$SBS18_cat, # Categorical
                      num.rep=1000)
# saveRDS(rp_SBS18, "results/random_forests/rp_SBS18")

rp_SBS40b <- rfPermute(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                       y = metabo_MAN$SBS1536A_cosmic, 
                       num.rep=1000)
# saveRDS(rp_SBS40b, "results/random_forests/rp_SBS40b")

rp_SBS22a_sro <- rfPermute(x = select(filter(metabo_MAN, country %in% c("Romania","Serbia")), c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                           y = filter(metabo_MAN, country %in% c("Romania","Serbia"))$SBS22_cat, # Categorical
                           num.rep=1000)
# saveRDS(rp_SBS22a_sro, "results/random_forests/rp_SBS22a_sro")

rp_SBS22b_sro <- rfPermute(x = select(filter(metabo_MAN, country %in% c("Romania","Serbia")), c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                           y = filter(metabo_MAN, country %in% c("Romania","Serbia"))$SBS1536I_cosmic_cat, # Categorical
                           num.rep=1000)
# saveRDS(rp_SBS22b_sro, "results/random_forests/rp_SBS22b_sro")

rp_DBS2 <- rfPermute(x = select(filter(metabo_MAN,!is.na(DBS2)), c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                       y = filter(metabo_MAN,!is.na(DBS2))$DBS2_cat, # Categorical
                       num.rep=1000)
# saveRDS(rp_DBS2, "results/random_forests/rp_DBS2")

rp_ID1 <- rfPermute(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                       y = metabo_MAN$ID1, 
                       num.rep=1000)
# saveRDS(rp_ID1, "results/random_forests/rp_ID1")

rp_ID5 <- rfPermute(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                       y = metabo_MAN$ID5, 
                       num.rep=1000)
# saveRDS(rp_ID5, "results/random_forests/rp_ID5")

rp_ID8 <- rfPermute(x = select(metabo_MAN, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                       y = metabo_MAN$ID8, 
                       num.rep=1000)
# saveRDS(rp_ID8, "results/random_forests/rp_ID8")

# Burdens
rp_SBSburden <- rfPermute(x = select(filter(metabo_MAN_burden_rf, !is.na(mutburden_sbs)), c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                          y = filter(metabo_MAN_burden_rf, !is.na(mutburden_sbs))$mutburden_sbs, 
                          num.rep=100)
# saveRDS(rp_SBSburden, "results/random_forests/rp_SBSburden")
rp_DBSburden <- rfPermute(x = select(filter(metabo_MAN_burden_rf, !is.na(mutburden_dbs)), c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                          y = filter(metabo_MAN_burden_rf, !is.na(mutburden_dbs))$mutburden_dbs, 
                          num.rep=100)
# saveRDS(rp_DBSburden, "results/random_forests/rp_DBSburden")
rp_IDburden <- rfPermute(x = select(filter(metabo_MAN_burden_rf, !is.na(mutburden_id)), c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                          y = filter(metabo_MAN_burden_rf, !is.na(mutburden_id))$mutburden_id, 
                          num.rep=100)
# saveRDS(rp_IDburden, "results/random_forests/rp_IDburden")

# SVs
rp_SVA <- rfPermute(x = select(metabo_MAN_SV, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                   y = metabo_MAN_SV$SBSSVA, 
                   num.rep=100)
# saveRDS(rp_SVA, "results/random_forests/rp_SVA")
rp_SVB <- rfPermute(x = select(metabo_MAN_SV, c(all_of(metabolites), "age_diag", "sex", "bmi","batch", "acq_order")), 
                    y = metabo_MAN_SV$SBSSVB_cat, # Categorical
                    num.rep=100)
# saveRDS(rp_SVB, "results/random_forests/rp_SVB")
